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Abstract 

Compressed Sensing (CS) has been applied in dynamic IVlagnetic Resonance Imaging (IVIRI) to accelerate the data 
acquisition without noticeably degrading the spatial-temporal resolution. A suitable sparsity basis is one of the key 
components to successful CS applications. Conventionally, a multidimensional dataset in dynamic MRI is treated as a series 
of two-dimensional matrices, and then various matrix/vector transforms are used to explore the image sparsity. Traditional 
methods typically sparsify the spatial and temporal information independently. In this work, we propose a novel concept of 
tensor sparsity for the application of CS in dynamic IVIRI, and present the Higher-order Singular Value Decomposition 
(HOSVD) as a practical example. Applications presented in the three- and four-dimensional IVIRI data demonstrate that 
HOSVD simultaneously exploited the correlations within spatial and temporal dimensions. Validations based on cardiac 
datasets indicate that the proposed method achieved comparable reconstruction accuracy with the low-rank matrix 
recovery methods and, outperformed the conventional sparse recovery methods. 
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introduction 

Dynamic MRI (dMRI) plays a vital role in many clinical 
applications, such as cardiac, perfusion and functional brain 
imaging. In these applications, high spatial-temporal resolution is 
desired to reveal anatomical details and physiological dynamics. 
Conventionally, the data is acquired in chronological order 
adhering to Nyquist sampling theorem, making MRI a relatively 
slow imaging modality. Routine methods speed up the MRI 
acquisition using a combination of fast gradient and Radio 
Frequency (RF) pulsing with full A-space sampling [1,2]. However, 
owing to hardware and physiological constraints, achieving high 
spatiotemporal resolutions with hardware intensive sequences is 
technologically challenging. 

Instead of increasing the data sampling rate, various approach- 
es, including Compressed Sensing (CS) [3], have attempted to 
reconstruct full field-of-view (FOV) images from sub-Nyquist 
acquisitions. CS has been recently apphed to MRI to accelerate 
the data collecting process. The pioneering work of applying CS to 
MRI to accelerate the data collecting process can be found in 
[4,5]. CS states that a faithful reconstruction of the signal is 
achievable with a sampling rate far lower than the Nyquist limit, 
provided that the signal has a sparse representation in some 
transform basis (called the 'sparsity basis'), which must be 
incoherent with the sensing matrix (i.e., Fourier transform in 
MRI) [3,6]. In static MRI and dMRI, the incoherence between 
the sensing basis and the sparsity basis can be achieved by 
randomly acquiring data in the A:-space or k-t space [3,7]. Both the 
predefined sparsity bases [8] and the data-dependent (also called 
data-derived) transforms [9,10] have provided successful recon- 
structions in static MRI applications. 



CS has also been applied to dMRI, where the data sets are 
naturally higher-order tensors (for instance, a third-order tensor for 
a cine MRI and a fourth-order tensor for a volume dMRI). 
Conventionally, 2D/1D sparsity bases were used to account for 
the spatial and temporal sparsity. When the method k-t SPARSE 
[4] was applied to the cine cardiac data, the 2D wavelet transform 
was first applied in the spatial domain, followed by the 1 D Fourier 
transform along the temporal dimension. The non-linear conju- 
gate gradient algorithm [1 1] was then used to reconstruct the 
sparsity coefficients. This is a practical and straightforward 
extension of the SPARSE MRI [8] as used in the static scenario. 
However, using 2D wavelet transforms may generate smooth/ 
blurry reconstructions at the image boundaries. Alternatively, the 
k-t FOCUSS method [12,13] applied different transforms to 
sparsily diverse MRI signals and explored the temporal sparsity by 
employing Principle Component Analysis and Fourier transform 
for the aperiodic and periodic/pseudo-periodic data, respectively. 
Then the recursively weighted minimum norm reconstruction 
algorithm (called TOCUSS') [14,15] was used to reconstruct the 
sparsity coefficients. Also using the FOCUSS algorithm, the k-t 
ISD [16] improved the CS reconstruction by exploiting the 
support information from the x-f space. Recent methods studied 
the anatomical structures or features [17,18] to further improve 
the reconstruction. Extending the application of sparsity, theoret- 
ical works [19-25] have investigated the low-rank matrix 
completion/recovery for more efficient signal recovery. The 
applications of the low-rank matrix structure have demonstrated 
merits in exploring the spatial-temporal signal redundancy in 
dMRI. For example, the methods described in [26-34] used sparse 
samphng schemes for data acquisition, and then generated basis 
functions for low-rank regularisation or to model the dMRI 
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signals. The function bases in methods [26-28] were tailored from 
the training data of the objects, they were more capable of 
capturing the correlations among the dynamic image series. The 
quality of the reconstructions achieved by these methods, however, 
relied heavily on the quality of the training data. Some other 
methods [29-34] used the combination of sparse sampling and 
low-rank regularisation without training data. 

Essentially, most of the existing CS-dMRI methods intend to 
use 2D/ID transforms to solve 3D or even higher-dimension 
problems. They either treat the 3D/4D data as a series of 2D 
images and then employ 2D/1D sparsilying transforms to explore 
spatial/ temporal sparsity [4,12,13] or, unfold the higher-order 
dataset into a 2D matrix to explore the spatiotemporal redun- 
dancy [29-31,35]. Intuitively, using matrix/vector transforms in 
dMRI data, being a higher-order tensor in nature, may not 
simultaneously explore the inherent data redundancy. To inves- 
tigate the possibilities of preserving the higher-dimensional data 
format, this work proposes a novel concept of tensor sparsity for 
dMRI. Inspired by a recent application of the second-order Singular 
Value Decomposition (SVD) [9, 1 0] in exploiting in-plane sparsity, 
the Tucker model based Higher-order Singular Value Decompo- 
sition (HOSVD) [36,37] was employed as a practical example for 
the current investigation. Tensor sparsity or tensor rank, is a 
powerful multidimensional signal processing tool that has been 
successfully applied in various areas. For instance in the area of 
pattern recognition/ computer vision, HOSVD has been used to 
extract the features of the training dataset to recognise/ classify 
future images (such as face verification) [38,39]. Recently, a low n- 
rank tensor approach has also been successfully applied to dMRI 
to achieve high quality image reconstruction for parallel and 
dMRI [33]. Instead of regularising the global low-rank structure, 
improved reconstruction accuracy and resolution were achieved 
by exploiting the local low-rank structure for multidimc'nsioiml 
MR signals, where the unknown values of the image matrices were 
locally estimated by considering the correlation among neighbour 
pixels or voxels [32,34]. Comprehensive reviews of the applica- 
tions of tensor decomposition, are provided in [40,41]. The 
HOSVD in the current study takes advantage of the fact that the 
signals in dMRI scenario are higher-order tensors. The presented 
approach sparsifies the dMRI signals in their original tensor 
format instead of the matrix format. Three experiments were 
designed to present the comparisons of the performances between 
this tensor sparsity basis and matrix transforms. In the first and the 
second experiments, the third-order SVD (3D-SVD) was used to 
sparsify the cine cardiac data (two spatial dimensions plus one 
temporal dimension). These experiments aim to compare the 
performance of the proposed method for pseudo-periodical data 
with two existing methods in dMRI. In the third experiment, the 
fourth-order SVD (4D-SVD) was applied as sparsity basis for the 
dynamic volume cardiac data (three spatial dimensions plus one 
temporal dimension), where the feasibility of the proposed sparsity 
basis in 4D application is demonstrated. 

The remainder of this article is organised as follows. Section 2 
explains the theoretical background of the proposed method. 
Section 3 describes the materials and methods used for validations. 
Section 4 presents the comparisons of the reconstruction results 
between the proposed method and the existing methods. Section 5 
discusses additional properties of the proposed sparsity basis. 
Section 6 concludes the contribution of this work. 

Theory 

In this section the general formulation of dMRI reconstruction 
in CS framework is first introduced. Then, the construction of a 



key component, the spasifying transform using tensor decompo- 
sition, is described. 

2.1. Formulation of Compressed Sensing in Dynamic iVlRI 
(CS-dMRI) 

To assist the discussion, the notations of scalars, vectors, 
matrices {.mond-order tensors) and tensors are denoted by 
lowercase letters (a, b, . . .), capital letters (A, B, ...) and caUigraphic 
letters {A, B, ...), respectively. Letter i and J are used to index row 
and column vectors, respectively. (A)j = Aj=aj, for example, 
denotes the jth column vector of matrix A. Hence, A= (Aj, A2, 
. . ./4j), where / is reserved for the index upper bounds, as is /. {A),j, 
also symbolised as fl^, denotes the element with a row index i and a 
column index J. 

Suppose an jWA-order tensor Ae C^' ^^^^ "' ^ is used to 
represent the spatial-temporal behaviour of the imaged object. 
Without losing generalit)', the first M = ,N-l dimensions of the 
tensor are used to describe the spatial information (for example 
M= 2 for 2D slice or M= 3 for 3D volume), which is collected at 
/^r time instances. The CS-dMRI problem can be solved using the 
following optimisation procedure: 

Minimise : \\H'iA)\\n 

(1) 

s.t.\\a>p(A)-y\\<e 

where y is the k space measurements collected from the MRI 
scanner; s represents the data-fidelity tolerance between the 
optimisation result and the measurements; *P is a transform that 
sparsifies the tensor A (the imaging object), and Oj- is a 
combination of operations, that is, the 2D Fourier transform for 
the in-plane data followed by a random under-sampling. 

Equation (1) minimises the /o-norm to enforce the sparsity of the 
object A, and uses the /2-norm as a constraint to guarantee the 
data-fidelity in the sampling domain. The optimisation problem in 
equation (1) is NP-hard (Non-deterministic Polynomial-time hard). 
The common solution for this problem is to relax the Iq norm to li 
norm, its nearest convex constraint [42]. Thus, the problem in 
equation (1) can be restated as: 

Minimise : \\^(A)\\i 
s.t.\\^AA)-y\\<e 

However, as has been extensively studied [43-48] , replacing the 
norm with an tp quasi-norm (()</;< 1) problem can reduce the 
amount of measurements needed for reconstruction or, can 
improve the reconstruction quality given the same amount of 
measurements. Therefore, in this work, we adopt the Ip norm as a 
constraint to enforce the sparsity of the images. Thus, the NP-hard 
problem in equation (1) can be replaced by solving a problem as 
follows: 

Minimise : \\'i'(A)\\': 

(3) 

s.t.\\Q>F(A)-y\\<e 

where 0<p<l. Section 3 will describe in detail the algorithm 
adopted to solve the non-convex problem in the form of equation 
(3). 
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2.2. Construction of the Sparsity Basis W using Higher- 
order Singular Value Decomposition 

In this section, the general framework of HOSVD [37] and the 
applications of HOSVD as sparsifying transform in CS-dMRI will 
be introduced. Several higher-order tensor operations wiU be 
introduced first to pave the way for the discussion of HOSVD. 
The HOSVD sparsity basis was obtained from the inverse Fourier 
transform of the zero-fdled under-sampled k space (denoted as .4o), 
therefore no training data is required in this method. 

Definition 1: Matrix unfolding-unfolding a tensor into matrices 
[37]. 



For the JWA-order tensor Aoe C' 



/| X/2 X ... X/^r 



the unfolded matrix 



Ao(„) contains the element ao,|,2 ,„ at the position with a row 
number i„ and a column number equal to 

(in+\ — ^)In + 2ln + 'i---lNhh---In-\ + 

('« + 2 — ^)In + iIn + A---Ilh---In-\ + ■■• + 

( //V - 1 )/l /2 . . . /„ _ 1 + ( /2 - 1 )/3 /4 • • • ^ - 1 + ■ ■ ■ + 4 - 1 



Figure 1 (a) exemplifies the process of unfolding a 3D tensor. There 
are three matrix representations (horizontal, lateral, and frontal) of 
the 3D tensor A^e C^' '''^'''^ in which all the slices are stacked one 
after another. The lateral matrix representation ^o(i)£ C^'^'^^^ is 
defined as [^0(l)]{,|_i)/3 + ,3 ,2 =flo,|,j,j; the frontal matrix represen- 
tation ^0(2)2 C'^''''^' is defined as [^o(2)l(,j_i)/|+,-,,„ =ao,,,j,3; and 
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Figure 1. Visualisation of matrix unfolding of a tensor, (a) 

Visualises the matrix unfolding of a third-order tensor and, (b) visualises 
a third-order tensor multiplied by matrix. 
doi:1 0.1 371/journal.pone.0098441 .gOOl 



the horizontal matrix representation ^o(3)E C^'^^^^' is defined as 
[^o(3)l(„_i)/2 + ,2,„=ao,„,,,. 

Definition 2: Multiplication of a higher-order tensor by matrices 
[37]. 

The n-mode product of the tensor AqS C'^^'-^ '"^'^ by a matrix 
UeC'"^'", denoted as ^x„(7 is an (IjXl2X...xI„_]Xj„xI- 
„+jX... xlj^ tensor, of which the entries are given by 

^■^0 ^«^^'l'2-'«-l'«'«+l-W ~ X]%/2...'„_1'«'„+1. 



Figure 1(b) visualises the multiplication of a 3D tensor by 
matrix, where S = C x 1 Ki x 2 K2 x 3 F3 {B = C''' ■'\C = 
i^iixijxhy figure 1(b) we can see that there are three 
multiplications of a 3D tensor by matrix. The 1-mode product of 
the tensor C = C^' ^'^^^^ by a matrix V] is defined as 

(CxiKi ),,,,,, = I] c„ ,2,3 V/,,,, which is a J1XI2XI1, .sized tensor; 
'1 

the 2-mode product of C = C^' by a matrix V2eC'^^^'^ is 

clef 

defined as (C x 2 F2)„,-,,3 = J2 '^'ii'ih^jih, which is a Iixjjxl^ 

ii 

sized tensor; similarly, the 3-mode product of C = C^' ^'- '^'^ by a 

matrix V^e C''^ " ^' is defined as (C x 3 K3),-,,^3 Y, ^hhh ''hh, 

h 

which is an I\ x IjX sized tensor. 

With these two operations defined, any Mh-order tensor 



/, X /2 X ... X /at 



Aoe C 
framework, by 



can now be decomposed, in the HOSVD 



flo 



..«S?US,...«^^L (4) 



h h JN 



where M*"j^,« = 1,2,...,A^, are the entries of the unitary matrices 
Ui,,n = 1,2,. ..,N, and Sj^j^^,„jf^ are the entries of Se C^i xh -^'N 
which is a complex tensor of size Ii x I2X ... x If^ . To facilitate the 
understanding of the properties of HOSVD, we first retrospect to 
the matrix SVD, which was used as a sparsity basis for static MRI 
[9, 1 0] . For any complex matrix Me C^' ''^h^ (^^n decompose it 
into product as 



M=USV" = Sx^Ux2V" = SxxUx X2U2 



(5) 



where U„,n = \,2, are /„ x/„,w=l,2, sized unitary matrices, and 5 
is an I\ X I2 sized matrix with the properties of [37]: 

(i) pseudo-diagonality; S = diag(a\,a2,...,Omm{ii,h)) 
(u) ordering! (T^ > C72 ^ ... ^ (Tjuin(/i,/2) 

>0 

where u,- are called the singular values of M. 
Likewise, in higher-order situation, we can decompose any 
complex Mth-ovAtr tensor Aae C'^^'^^ '"^'^ as 



Ao = SXiUlX2U2X^...XfiUN, 



(6) 



where the unitary matrices U„e C'" ,n = \,2,...,N are called the 
n-mode singular matrices. Tensor ^^h^h'^ ■■■'^'n j^^j 
following properties: 
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(i) all-orthogonality: two sub-tensors >Si„=a and Si^^p are 
orthogonal for all possible n, a. and P subject to (X P, 
which means '(5,„=(n5,„=/() = 0 when a p, 

(ii) ordering: ||<Si„ = i|| > ||<Si„=2|| >-> ||>5i„=/J| >0 for aU possi- 
ble n, 

where the Frobenius-norms I|5,„=/||, symbolised by <7^"\ are called 
the n-mode singular values of Aq- 

As demonstrated in [37], given a Mh-oider tensor Afi, the n- 
mode singular matrix U„ in equation (6) is actually the left singular 
matrix of the correlated n-mode matrix unfolding of (as per 
Definition 1 and 2). Therefore the computation of the HOSVD in 
equation (6) eventually leads to jV different matrix SVD operations 
on the unfolded tensor. Therefore, the tensor S can be computed 
as 

5 = AXlC/,^X2f/f X3...Xjvi7^ (7) 

For example, Ui can be obtained by performing the matrix 
SVD on the 1-mode unfolding matrix ^o(l) as: 

Si = C/i^0(i)Fi (8) 

Generally, U„,n= l,2,...,N, can be obtained by performing the 
matrix SVD on the n-mode unfolding matrix ^o(n) as: 

S„=UnAo^n)Vn. (9) 

With the unitary matrices obtained, we can then construct the 
tensor sparsifying transform as: 

^>{Ao) = AoXiUfx2U^X3...X]^U§ (10) 

where the sparsity basis U„,n=\,2,...,N is obtained from the 
inverse Fourier transform of the zero-filled under-sampled k space 

The inverse sparsifying transform is then obtained as: 

'¥-\S)=SxiUiX2U2X3...Xf^Uf^ (11) 

Figure 2(a) visuahses the decomposition of a third-order tensor 



Ao=SxiUiX2U2X^U^ (12) 

The unitary matrices [/„,w= 1,2,3, in equation (12) can be 
obtained from equation (9). The properties of all-orthogonality 
and ordering [37] guarantee that most of the energy of tensor S 
accumulates around one vertex, and little energy distributes to the 
broad area away from this region. Therefore, tensor S has a 
sparse representation (refer to figure 2(a) for illustration). Likewise, 
figure 2(b) presents an example of the HOSVD in the fourth-order 
tensor case. It should be noted that the tensor S is shown in 
logarithmic scale to assist the presentation, because S is too sparse 
to be easily visible. It is clearly shown that in both 3D and 4D cases 
the coefficients with large values are highly concentrated in one 



voxel (light blue colour), while the vast majority of the elements in 
the S tensor are close to zero (deep blue colour). 

Materials and Methods 

To test the possibility of employing HOSVD as higher-order 
sparsifying transform in GS-dMRI applications, three experiments 
are designed: two cine cardiac MRI schemes and one dynamic 
volume cardiac MRI series. 

3.1. Datasets 

3.1.1. 3D-SVD; Application in cine cardiac MRI. Two sets 

of cine cardiac MRI data were used to validate the proposed 
method. The first dataset (Dataset A) was acquired at the 
University of Utah, which was used in the method k-t SLR [30] . 
70 frames of k-space were acquired on a 3T Siemens scanner with 
the spatial resolution of 90x190 (phase encoding xfrequency 
encoding). The cardiac data was obtained with a saturation 
recovery sequence (TR/TE = 2.5/1 ms, saturation recovery 
time= 100 ms). The second dataset (Dataset B) was acquired at 
Yonsei University Medical Center, which was used in the method 
k-t FOCUSS [12,13]. 25 frames of full k-space data was acquired 
using a 1.5T Philips system with an in-plane spatial resolution of 
256x256. The cine cardiac data was obtained using steady-state 
free precession (SSFP) sequence with a flip angle of 50 degree and 
TR = 3.45 msec. The FOV was 345 mmx270 mm. The slice 
thickness was 10 mm. A few frames from both Dataset A and 
Dataset B are shown in figure 3(a, b). 

3.1.2. 4D-SVD: Application in volume dynamic cardiac 
MRI. The third experiment investigated the possibility of 
employing the proposed HOSVD sparsifying transform for 4D 
dynamic cardiac MRI. The images of the this dataset (Dataset C) 
were of one subject, arbitrarily chosen from a total of 33 available 
subjects [49]. The measurements were acquired from a GE 
Genesis Signa MR scanner using the FIESTA protocol. The 
dimension of the subject data is 256x256x10x20 (phase 
encoding xfrequency encoding xz position X time). It is noted that 
Dataset C is in DIGOM (Digital Imaging and Communications in 
Medicine) format. Using the real-valued images, instead of the 
complex-valued k-space data, has made the reconstruction of the 
4D experiment easier. A few frames from Dataset C are shown in 
figure 3(c). 

3.2. Reconstruction 

3.2.1. Optimisation algorithm. The Ip quasi-norm in 
equation (3) poses a non-convex optimisation problem. Theoret- 
ical work [44,46] has demonstrated that this non-convex problem 
is solvable and, the local minima can be avoided in practice 
[43,50]. The applications in the medical imaging context 
[30,45,47,48,51-55], have already demonstrated the practicability 
and advances of non-convex optimisation. In this work, we adopt 
the algorithms used in [30,47] to solve the optimisation problem 
stated in equation (3). In [47], Chartrand used both wavelet 
transform and discrete gradient to enforce the sparsity of the MR 
images. In [30], Lingala et al. used the combination of rank 
property and signal sparsity for reconstruction. In this work we use 
only the HOSVD for sparsity enforcement. Therefore we herein 
briefly state the modified optimisation process as follows. 

We begin with the definition of a variable splitting operator: 

®iA) = Minimise : ||»P(B)||^-I-;6||.A-B||2 (13) 
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Figure 2. Visualisation of HOSVD. (a) Shows HOSVD on a third-order tensor and, (b) shows HOSVD on a fourth-order tensor. The left of (b) shows 
a four dimensional cardiac dataset denoted as A. The four dimensions are labelled as /,, Z^, is and /V The right of (b) presents tensor S and the unitary 
matrices U, U2 U3 and U4, that were obtained by performing HOSVD operation on A. S is also a fourth-order tensor, the dimensions of which are 
marked as /,, i2, 13 and Z^. 
doi:1 0.1 371 /journal.pone.0098441 .g002 



where /?>0 is a constant. It is noted that &{B) is forced to 
approach when /?^oo. 

We rewrite the problem in equation (3) as into its Lagrange's 
form as: 



Minimise : ||a)^.(^)-j;||2 + 



(14) 



where X is a constant to balance the weighting between the data 
fidelity and the signal sparsity. Then the splitting operator was 
applied on equation (14), arriving at: 



Minimise : \\Q>F(A)-y\\2 + P\\A-B\\ 



(17) 



To solve the sub-problem with respect to variable B, we fix 
variable A and apply ^-shrinkage operator to each pixel oi^(B). 
As explained in [47] the /i-shrinkage operator is executed as: 



S'^ih) = max[\b\ - a.^- ' ,Q}bl\b\ 



(18) 



Minimise : ||<l)f(^)-ji'||2 + ;.0(^) 
which can be expanded as: 



(15) 



Minimise : 



\\<fFiA)-y\\, + fi\\A-B\\, + MmB)\\'; (16) 



We can then solve the problem above by iteratively solving the 
variables A and B in turn. In this way the problem in equation 
(16) is decomposed into two simple sub-problems. The two sub- 
problems are decoupled, making it computational efficient. By 
setting fS—^co, the solution of equation (16) approaches that of 
equation (14). 

To solve the sub-problem with respect to variable A, we can fix 
variable B and adopt the conjugate gradient algorithm as used in 
[30]: 



To choose an appropriate value for the parameter [i, we 
initialised it with a relatively small value and then geometrically 
increased it as proposed in [56]. To enforce the data-fidelity, the 
residual of each sub-problem was added back to the data at each 
iteration as proposed in [57]. For a summary of the optimisation, 
please refer to figure 4. 

3.3. Comparison Validations 

The proposed method was compared with one of the recent 
low-rank image reconstruction methods, k-t SLR, and a classic CS 
method, k-t SPARSE. For fair comparison, we ensure that firstiy 
all the methods used the same sampling pattern of k-t space; 
secondly, the parameters for all the methods were adjusted 
appropriately so that both the signal to error ratio (SER) and the 
averaged signal intensity for all methods were optimised; and 
thirdly the optimisations for all the methods share the same 
stopping criterion, that is the optimisations ceased when the 
gradient magnitude of the object function reached 1x10 or the 
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' ' r -I /'• / / 



Figure 3. Several frames of the datasets obtained at different 
time instants (as indicated). From top to bottom; (a) Dataset A, (b) 
Dataset B and, (c) Dataset C. In (c), images obtained at four time instants 
(as indicated by fj) are presented row by row, and images obtained at 
fourz positions (as indicated by z, along axial direction (base-apex)) are 
presented column by column. The regions of interest of each dataset 
are marked within the red rectangles and, the regions of myocardial 
and blood pool used for averaged signal intensity comparison are 
marked with yellow and light blue colours, respectively. 
doi:1 0.1 371/journal.pone.0098441 .g003 

number of iterations reached 300. All the evaluations were 
implemented using Matlab 201 la (MathWorks, Natick, MA) on a 
Mac OS X Lion operation system, with a dual-core 2.4 GHz Intel 
processor and 4 GB of memory. The SER was calculated as: 

5£^=-101ogio^^7^^%^, (19) 

1 1 "^/"'Z I 

where A„s is the result of the reconstruction, Aj,M is the fuUy 
sampled dynamic images, and ||'||^ denotes the Frobenius norm. 
A greater SER value correlates to a better image quality. 

The method k-t SLR employs two regularisations: the low-rank 
structure and the sparsity of the signal. To exploit the low-rank 
structure, k-t SLR reshaped the 3D dataset into a large 2D matrix 
r. More specifically, the 2D images in a dynamic sequence were 
firstly vectorised and then concatenated to form the matrix F. To 
exploit the sparsity of the signal, the total variation (TV) was used 
as an extra regularisation. Moreover, instead of using convex 



Input: k-t space data^", parameters X, P; 

Initialisation: Ai = Of"'(v), b] = Q 

for n=\: the number of outer iterations do 

while the stopping criterion is NOT satisfied do 
m=\; 

update Am to Am+ 1 by solving equation ( 1 7); 

apply shrinkage operator 4'(Bm+i) = S''a('+'(Am+i) + b,,,); 

b„,^i = b„, + H^CBm+i); 

Am+l = Bni+i; 

m = m + 1 ; 

end 

P,,,; = P„ X INCREASE _PARAMETER: 

end 



Output: A 



Figure 4. Outline of the reconstruction algorithm. 

doi:10.1371/journal.pone.0098441.g004 

penalties to regularise the low rank property and the sparsity, k-t 
SLR adopted some of the recent algorithms on the non-convex 
regularisation [43,47,48] for the optimisation, further improving 
the reconstruction result. In [30], the combination of the 
constraints provided better image quality than the variants of 
the k-t SLR, which rely on either matrix SVD or TV constraint 
alone. Therefore in this work, we only compare the proposed 
method with k-t SLR, where both SVD and the TV regularisations 
were used in the optimisation. The method k-t SPARSE is a classic 
CS-dMRI method. It uses the wavelet transform (Daubechies 4 
was used as the mother wavelet in this work) for in-plane sparsity 
and the Fourier transform for temporal sparsity, assuming that the 
change of the heartbeat is periodical. AH the methods compared in 
this work are flexible to account for arbitrary non-Cartesian k- 
space sampling schemes; we adopt the radial trajectory with 
uniform angular spacing as used in [30]. The trajectory was 
randomly rotated with a small angle for each frame to implement 
random sampling. 

Results 

4.L 3D Application 

Figure 5 and 6 show the reconstruction of Dataset A at 
reduction factors 6 and 1 1 respectively. When the reduction factor 
was 6 (reduction factor n means only 1/n of the full A-space 
measurements were obtained), the SER values achieved by the 
proposed method, the k-t SLR and the k-t SPARSE, were 8.9 dB, 
8.7 dB and, 7.7 dB, respectively. Figure 5 shows the reconstruc- 
tion of Dataset A when the reduction factor was 6. As shown in 
figure 5(b), all the methods provided comparable averaged signal 
intensity for the blood pool area (normalised to the maximum grey 
level of the region of interest, figures 6, 7, 8, and 9 are normalised 
in the same fashion). However, when comparing at the myocardial 
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Figure 5. Reconstructions of Dataset A at reduction factor 6. (a) and (b) show the averaged normalised signal intensity at the myocardial and 
blood pool regions, respectively, and (c) shows the images (region of interest only) at the peak signal intensity of myocardial (the 54* frame, top row) 
and blood pool (the 14* frame bottom row). The left of (a) and (b) shows the averaged signal intensity of the fully sampled images (black line), k-t 
SLR reconstruction (blue line) and, the reconstruction of the proposed method (red line); the right of (a) and (b) shows the averaged signal intensity 
of the fully sampled images (black line), the k-t SPARSE reconstruction (blue line) and, the reconstruction of the proposed method (red line). 
doi:1 0.1 371 /journal.pone.0098441 .g005 



signal intensity, both the proposed method and k-t SLR obviously 
outperformed k-t SPARSE, especially at the frames where the 
averaged signal intensity changed rapidly (as marked with red 
arrows in figure 5(a)). The region of interest (as marked in 



figure 3(a)) of the 54* and the 1 4* frames, where the myocardial 
and the blood pool signal intensities reached their peak values, are 
presented on the top and the bottom rows of figure 5(c), 
respectively. In figure 5(c), it appears that Dataset A contains 
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Figure 6. Reconstructions of Dataset A at reduction factor 11. (a) and (b) show the averaged normalised signal intensity at the myocardial 
and blood pool regions, respectively, and (c) shows the images (region of interest only) at the peak signal intensity of myocardial (the 54'*^ frame, top 
row) and blood pool (the 14"^ frame bottom row). The left of (a) and (b) shows the averaged signal intensity of the fully sampled images (black line), 
k-t SLR reconstruction (blue line) and, the reconstruction of the proposed method (red line); the right of (a) and (b) shows the averaged signal 
intensity of the fully sampled images (black line), the k-t SPARSE reconstruction (blue line) and, the reconstruction of the proposed method (red line). 
doi:1 0.1 371 /journal.pone.0098441 .g006 



visible white noise, and some of the residual noise was maintained 
in the result of the method k-t SPARSE. The images recovered by 
the proposed method and the k-t SLR successfully supressed the 
white noise. Both the proposed method and the k-t SLR provided 



comparable overall image quality at the low reduction factor. 
When the reduction factor was 1 1 , the SER values achieved by the 
proposed method, the k-t SLR and the k-t SPARSE, were 8.0 dB, 
7.8 dB and, 6.4 dB, respectively. As shown in figure 6(b), all the 
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Figure 7. Reconstructions of Dataset B at reduction factor 6. (a) and (b) show the averaged normalised signal intensity at the myocardial and 
blood pool regions, respectively, and (c) shows the images (region of interest only) at the peak signal intensity of myocardial (the 1 3* frame, top row) 
and blood pool (the 24* frame bottom row). The left of (a) and (b) shows the averaged signal intensity of the fully sampled images (black line), k-t 
SLR reconstruction (blue line) and, the reconstruction of the proposed method (red line); the right of (a) and (b) shows averaged the signal intensity 
of the fully sampled images (black line), the k-t SPARSE reconstruction (blue line) and, the reconstruction of the proposed method (red line). 
doi:1 0.1 371 /journal.pone.0098441 .g007 



methods recovered comparable averaged signal intensity of the 
blood pool area. However, when comparing the myocardial area, 
the proposed method and the k-t SLR outperformed the k-t 
SPARSE more obviously, especially at the frames where the signal 
changes quickly (as indicated by the red arrows in figure 6(a)). The 



region of interest of the 54 and the 14'^ images are presented on 
the top and the bottom row of figure 6(c), respectively. The k-t 
SPARSE was severely affected by the white noise at this high 
reduction factor, while both the k-t SLR and the proposed method 
were still robust to noise. In k-t SLR, TV regularisation provided 
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Figure 8. Reconstructions of Dataset B at reduction factor 1 1 . (a) and (b) show the averaged normalised signal intensity at the myocardial and 
blood pool regions, respectively, and (c) shows the images (region of interest only) at the peak signal intensity of myocardial (the 1 3* frame, top row) 
and blood pool (the 24* frame bottom row). The left of (a) and (b) shows the averaged signal intensity of the fully sampled images (black line), k-t 
SLR reconstruction (blue line) and, the reconstruction of the proposed method (red line); the right of (a) and (b) shows the averaged signal intensity 
of the fully sampled images (black line), the reconstruction of k-t SPARSE (blue line) and, the reconstruction of the proposed method (red line). 
doi:10.1371/journal.pone.0098441.g008 



slightly better reconstruction for large contours or boundaries of 
the images. However, it also generated a cartoon-like /over- 
smooth effect on local fine details (also observed in [30]). This 
effect is more obvious at reduction factor 1 1 (see figure 6(c)). 
Compared with TV regularised k-t SLR, the proposed method 



provided sUghtly better reconstruction of local fine details (see the 
red arrows in figure 6(c)). The evaluation of all the methods based 
on Dataset A indicates that the proposed tensor sparsity basis 
outperformed the conventional matrix sparsity basis. Moreover, 
even when comparing with k-t SLR that combines the low-rank 
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Figure 9. The averaged normalised signal intensity achieved 
by the proposed method at reduction factor 11. (a) The 

myocardial signal intensity of the fully sampled images and the 
reconstructed images provided by the proposed method; (b) the blood 
pool signal intensity of the fully sampled images and the reconstructed 
images provided by the proposed method. 
doi:1 0.1 371/journal.pone.0098441 .g009 

matrix recovery and the sparsity constraint, the proposed method 
was still able to provide comparable overall reconstruction 
accuracy. 

Figure 7 and 8 show the reconstruction of Dataset B provided 
by all the methods at reduction factors of 6 and 1 1, respectively. 
When the reduction factor was 6, the proposed method, the k-t 
SLR and the k^t SPARSE achieved the SER values of 12dB, lOdB 
and, 9.9dB, respectively. The averaged signal intensity compar- 
ison, as shown in figure 7(a, b), demonstrates that the proposed 
method was more capable of capturing the dynamic features of the 
signal (see the red arrows in figure 7(a, b)) than k-t SPARSE. The 
13* and the 24* frames (region of interest only, as marked in 
figure 3(b)), where the peak averaged signal intensity of myocardial 
and blood pool areas were reached, are presented on the top and 
the bottom rows of figure 7(c), respectively. As shown in figure 7 
(c), Dataset B contains more local details than Dataset A and, it 
has litde visible white noise. All the methods succeeded in 
recovering the coarse features of Dataset B; meanwhile, the 
proposed method and the k-t SLR captured more fine details (see 
the red arrows in figure 7(c)). When the reduction factor was 1 1, 
the SER values achieved by the proposed method, the k-t SLR and 
the k-t SPARSE, were 10.4 dB, 9.5 dB and, 8.4 dB, respectively. 
The averaged signal intensity of the myocardial and the blood pool 
was compared in figure 8(a, b). The proposed method achieved 
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Figure 10. The reconstruction of Dataset C (region of interest 
only) achieved by the proposed method at reduction factor 1 1 . 

(a) Presents several fully sampled images, and (b) presents the 
corresponding reconstructed images. In (a), images obtained at four 
time instants (indicated by f) are present row by row; and images 
obtained at four z positions (indicated by z along axial direction (base- 
apex)) are presents column by column. Likewise, (b) presents the 
reconstructed images at the corresponding time instants (indicated by 
f) and z positions (indicated by z). 
doi:10.1371/journal.pone.0098441.g010 

comparable reconstruction with the k-t SLR and, better overall 
reconstructions as compared to k-t SPARSE. And the visual 
evaluation in figure 8(c) shows consistent results with those of the 
averaged signal intensity comparison, as indicated by the red 
arrows. The quantitative and visual evaluations of Dataset B were 
also consistent with those of Dataset A. 

4.2. 4D Application 

As shown in figure 2, the HOSVD method can be applied 
straightforwardly to higher order datasets. In this work, we present 
the application of HOSVD in the dynamic volume cardiac 
imaging, where the dataset is a 4D tensor. At reduction factor 1 1 , 
the SER of the reconstructed 4D images achieved by the proposed 
method was 12.1 dB. The averaged signal intensity at the 
myocardial and the blood pool areas is presented in figure 9. As 
illustrated in figure 9, at the high reduction factor of 1 1 the 
proposed method was still able to recover the dynamic features of 
the signal without noticeable error. Several fuUy sampled and the 
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reconstructed images (region of interest only, as marked in 
figure 3(c)) are shown in figure 10(a) and (b), respectively. As 
shown in figure 10(b), the proposed method successfully recovered 
the coarse features of the object and, most of the fine details were 
also recovered, which demonstrated the feasibility of the proposed 
sparsifying transform for the 4D application. 

Discussions 

5.1. The Tucker Model Based HOSVD 

This work takes the Tucker model based HOSVD as an 
example to demonstrate the potential of tensor decomposition in 
the exploration of higher-order signal sparsity. The Tucker model 
based HOSVD decomposes a dense tensor into a sparse tensor 
multiplied by matrices along individual modes (as shown in 
figure (1-2)). The k-t SLR actually used solely the mode-2 unfold of 
the tensor structure to explore the low rank properties. However, 
this work does not explore the low-rank structure of the reshaped 
tensor. Instead, it explores the sparsity in a tensor structure. In 
addition to HOSVD, th(;rc are a broad range of tensor 
decomposition techniques for future investigation, such as the 
CANDECOMP/PARAFAC decomposition [58,59] and its vari- 
ants, which can be used to explore the tensor rank minimisation. 

5.2. Computational Cost 

In this work, when the same stopping criteria was set, the 
computation time for the proposed method in the 3D application 
was, on average, 28 and 29 minutes for Dataset A and B 
respectively. The k-t SLR method used 2 1 minutes for Dataset A 
and 29 minutes for Dataset B. As for the k-t SPARSE, it took 
29 minutes on average for Dataset A and 40 minutes for Dataset 
B. In dealing with third-order tensor, the proposed method 
performs SVD three times (once per equation (8-10)), while k-t 
SLR needs only one SVD computation. The proposed method 
involves only one regularisation, while the k-t SLR involves two 
regularisations. Therefore, though the computing time of the 
HOSVD basis function is three times that of the SVD basis 
function, the overall optimisation time of the proposed method 
was only approximately 30% more than that of the k-t SLR. The 
sparsifying transform in the k-t SPARSE involves multiple times of 
wavelet transforms for each frame and one Fourier transform for 
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Conclusion 

This work proposes a novel concept of tensor sparsity for 
Compressed Sensing in dynamic MRI, and presents the Tucker 
model based Higher-order Singular Value Decomposition as a 
practical example. Tlu; tensor d(;composition based method 
derives the sparsity basis adaptively and directly from the zero- 
filled under-sampled k-t space measurements, and does not require 
extra scan time to obtain training data. The proposed tensor 
sparsity basis provides improved image reconstruction quality 
when compared to the classic sparsity basis. The reconstruction 
quality is similar to that with a stronger constraint-low rank 
property of matrix. 
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